function [irfa] = mk_irfa(Betat,Cm,nh)
% PURPOSE: generates IRF array (nn x nn x (nh+1)) 
% -----------------------------------------------------
% USAGE: [irfa] = mk_irfa(Betat,Cm,nh)
% where: 
% Betat = (N x N*nlags) VAR coefficient matrix=[A1 A2 ... Ap]
% Cm    = (N x N) Identification Matrix (e.g. Cholesky Matrix)
% nh = Horizon for IRFs
% NOTE: no dimension checks imposed
% -----------------------------------------------------
% MODEL : y(t) = A1*y(t-1) +...+ Ap*y(t-p) + Cm*e(t)
% e(t) = structural error
% -----------------------------------------------------
% RETURNS: 
% irfa = (N x N x (nh+1)) 3Dim Array with IRF coefficients
% nkk= N*nlags+nd
% The dimensions of irfa are then
% (number of outputs) � (number of inputs) � (length of t)
% and irfa(:,j,:) gives the response to an impulse disturbance entering the jth input channel. 
% -----------------------------------------------------
% written by:
% Andrea Tamoni
% andrea.tamoni@phd.unibocconi.it
% � Dec 2007

N=size(Betat,1);
nlags=size(Betat,2)/N;

irfa=zeros(N,N,nh);
%SET IMPULSE RESPONSES : j = 1
irfa(:,:,1)=Cm;

%Companion
Mm = var_compn(Betat);

Mm1=eye(N*nlags);
%IMPULSE RESPONSES : j > 1
for i=1:nh-1
    Mm1=Mm1*Mm;
    irfa(:,:,i+1)=Mm1(1:N,1:N)*Cm;
end